
clear
clear all
clear mata
clear matrix
set more off
set matsize 11000
set maxvar 30000
cap log off
capture log close
set emptycells drop
pause on


***********************************************
* USER DEFINE FILEPATH FOR REPLICATION FOLDER *
***********************************************

global path "C:\Users\wb520443\Dropbox\Research Projects\Global pollution\2_analysis\Replication"



***********************************************


global opts		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))  ) star(* .10 ** .05 *** .01) noobs nocons

global opts_fs		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))) star(* .10 ** .05 *** .01) noobs nocons


global input "$path\data"
global figdat "$path\figdat"
global figures "$path\figures"
global tables "$path\tables"



cap mat define countryfigure=J(105,10,.)
cap mat define countryfigure_median=J(105,10,.)
cap mat define countryfigure_max=J(105,10,.)
		local row=1
		
	import delimited "$input/dhs_cluster_wealth_pm_anthro.csv"

	
	g west_africa=1 if country=="Senegal"|country=="Sierra Leone"|country=="Nigeria"|country=="Morocco"|country=="Mali"|country=="Liberia"|country=="Ghana"|country=="Guinea"|country=="Cote d'Ivoire"|country=="Chad"|country=="Cameroon"|country=="Burkina Faso"|country=="Benin"|country=="Togo"|country=="Namibia"|country=="Gabon"
	
	destring pm_mean, replace force
	
	rename country country2
	destring pm_anthro_mean, replace force
	cap encode country, g(country_code)
		label save using "$input/labels.do", replace
	merge m:1 country2 using "$input/country_population.dta"
	
	replace country_population=1.601e+09 if country2=="Pakistan"
	replace country_population=26378275 if country2=="Cote d'Ivoire"
	replace country_population=89561404 if country2=="Democratic Republic of the Congo"
	replace country_population=6591600 if country2=="Kyrgyz Republic"
	replace country_population=59734213 if country2=="Tanzania"
	replace country_population=1318442 if country2=="Timor"
	
	
	replace country_population=1
	
	encode svy_id, g(id)
	

	//Mean, mean_ss, and max refer to how PM25 across years and months is collapse (e.g. max is the max pollution level for a RWI point across all years and months in the sample).
	forvalues i=1(1)103{
		
		cap reghdfe pm_mean wealth_n  if country_code==`i', absorb(i.id) 
			cap mat def countryfigure[`row',1]=_b[wealth_n]
			cap mat def countryfigure[`row',2]=_se[wealth_n]
		cap reghdfe pm_anthro_mean wealth_n if country_code==`i', absorb(i.id) 
			cap mat def countryfigure[`row',3]=_b[wealth_n]
			cap mat def countryfigure[`row',4]=_se[wealth_n]
		
		local row=`row'+1
	}	
	

	
clear
			
		svmat2 countryfigure

		drop if countryfigure1==.
		
		g country_code=_n
			run "$input/labels.do"
		
		g HI95_mean=(1.96*countryfigure2)+countryfigure1
		g LI95_mean=countryfigure1-(1.96*countryfigure2)

		g HI99_mean=(2.06*countryfigure2)+countryfigure1
		g LI99_mean=countryfigure1-(2.06*countryfigure2)
		
		
		g HI95_mean_ss=(1.96*countryfigure4)+countryfigure3
		g LI95_mean_ss=countryfigure3-(1.96*countryfigure4)

		g HI99_mean_ss=(2.06*countryfigure4)+countryfigure3
		g LI99_mean_ss=countryfigure3-(2.06*countryfigure4)
		
				
		set scheme plotplain
		
		
		rename (countryfigure1 countryfigure3 countryfigure5  ) (mean_beta mean_ss_beta max_beta )
		
		sort mean_beta
		g n=_n
		g n2=n+0.1
		
		

		label val country_code country_code
			decode country_code, g(country)
			replace country="IND/PAK" if country=="IND_PAK"
			labmask n, values(country)
			
		
		

		
		twoway 	(rspike HI99_mean LI99_mean n , lc(gs10) horizontal) ///
					(scatter n mean_beta  if mean_beta<-1.9, mc(cranberry%60) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-2  & mean_beta<-0.9, mc(cranberry%60) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-1 &  mean_beta<0, mc(cranberry%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>0  & mean_beta<1, mc(navy%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>1  & mean_beta<2, mc(navy%60) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>2  , mc(navy%90) ms(D) msize(large) ///
					yti("") ysc(noline axis(1)) ylabel(,nolab alternate valuelabel angle(0) labsize(large) nogrid noticks) ysc(noline)  xlabel(,nogrid labsize(large)) xtitle("Correlation between wealth" "and average pollution",size(vlarge)) ysize(100) xsize(30) legend(off) xline(0))
					
					
					graph export "$figures/country_figure_mean_collapsed_DHS.png", as(png) replace
			
		

clear
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	